//Course Project
//Olekhno N. 2103/2
//12.12.12

#include <iostream>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <math.h>
#include <graphics.h>

const int NS=100000;                  //number of steps in modeling
const int N=60;                       //number of layers in body

using namespace std;

main()
{
  double Fx[N], Fy[N], dFx[N], dFy[N]; //forces
  double x[N], y[N], dx[N], dy[N];     //coordinates
  double kx=1, ky=4;                   //material properties; actually here k=k/m
                                       //where m is a mass of layer
  float dt=0.1;                        //step-time of modeling
  int a=10;                            //distance between layers
  int i,j,k;                           //countes
  int width, height;                   //window parameters
  char q;                              //used in user dialog                              
  
  y[0]=0;                              //initialization
  x[0]=0;
  for (i=1; i<N-1; i++)
  {
     y[i]=0;
     Fy[i]=0;
  }
  for (i=1; i<N-1; i++)
  {
     x[i]=0;
     Fx[i]=0;
  }
  
  Fx[0]=50;
  Fy[0]=20;
  Fx[30]=-20;
  Fy[30]=15;
  width=1000;
  height=700;
  
  cout<<"Would you like to change settings or start with defaults? y/n"<<endl;
  cin>>q;
  if (q=='y' or q=='Y')
  {
     cout<<"Enter timestep of modeling dt, sec <0.1>"<<endl;
     cin>>dt;     
     cout<<"Enter distance between layers <10>"<<endl;
     cin>>a;     
     cout<<"Enter kx <1>"<<endl;
     cin>>kx;
     cout<<"Enter ky <4>"<<endl;
     cin>>ky;
     cout<<"Enter the value of horizontal component of force on zero layer <50>"<<endl;
     cin>>Fx[0];
     cout<<"Enter the value of vertical component of force on zero layer <20>"<<endl;
     cin>>Fy[0];
     cout<<"Enter the value of horizontal component of force on thirtieth layer <-20>"<<endl;
     cin>>Fx[30];
     cout<<"Enter the value of vertical component of force on thirtieth layer <15>"<<endl;
     cin>>Fy[30];
  }
  
  initwindow(width,height, "Oscillators");
  setcolor(RED);
  
  for (k=0; k<NS; k++)                     //step of modeling
  {
     for (j=0; j<300; j++)
     {
        for (i=0; i<N-2; i++)              //calculating new coordinates
        {
           dx[i]= 0.5*Fx[i]*dt*dt;
           x[i]=x[i]+dx[i];
           dy[i]= 0.5*Fy[i]*dt*dt;
           y[i]=y[i]+dy[i];
        }
       // x[0]=x[1];   //+in previous cycle i=1 is equal to a case with two
       // y[0]=y[1];   //free borders
        dx[N-1]=0;
        dy[N-1]=0;
     
        Fx[0]=-1*kx*(2*x[0]-x[1]);      //calculating new forces
        Fy[0]=-1*ky*(2*y[0]-y[1]);
        for (i=1; i<N-3; i++)
        {
           Fx[i]= -1*kx*(x[i]-x[i-1])-kx*(x[i]-x[i+1]);
           Fy[i]= -1*ky*(y[i]-y[i-1])-ky*(y[i]-y[i+1]);
        }
        Fy[N-2]=-1*ky*(y[N-2]-y[N-3]+y[N-2]);
        Fy[N-1]=0;
        Fx[N-1]=0;
     }
                      
     {                                  //visualization of movement
        Sleep(40);                      //25 frames per sec
        setcolor(WHITE);
        bar(0,0,1000,700);
        setcolor(RED);
        for (i=0; i<N-1; i++)
        {
           moveto(100+10000*x[i], a*(i+1)+1000*y[i]);
           lineto(100+10000*x[i]+200, a*(i+1)+1000*y[i]);
        }
     } 
  }
  
  getch();
  closegraph();
  system ("pause");
  return 0;                               
}
